#!/usr/bin/perl
# 2012/12/12 Calculate chromosome combined distance curve and output fitting curve paramter file
# 2012/12/10 Calculate raw data's distance curve

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 8 or $ARGV[0] eq '--help'){
	die "Usage : $0  -i [association map file] -o [output average file] -p [output parameter file] -g [output graph pdf name (don't use %)]\n";
}

my %opt;
getopts("i:t:r:o:p:g:", \%opt);
my $FILE_map = $opt{i};
my $FILE_average = $opt{o};
my $FILE_parameter = $opt{p};
my $FILE_graph = $opt{g};


#---------------------------------------
# read association amp
#---------------------------------------
my %SUM;
my %NUM;
my ($SUM_chr1_chr2, $SUM_chr1_chr3, $SUM_chr2_chr3) = (0,0,0);
my ($INTER_chr1_chr2, $INTER_chr1_chr3, $INTER_chr2_chr3) = (0,0,0);
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";

my $flagOne = 0;
my $Resolution = 0;
my %MAX_distance = ('chr1' => 0, 'chr2' => 0, 'chr3' => 0);

while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id1, $id2, $value) = split /\t/;
	my ($chr1, $pos1, $chr2, $pos2);
	if($id1 =~ m/(chr\d):(\d+)\-(\d+)/){
		$chr1 = $1;
		my $start = $2;
		my $end = $3;

		if($flagOne == 0){
			$Resolution = $end - $start + 1;
			$flagOne = 1;
		}

		$pos1 = int(($start + $end)/2);
	}else{
		die "cannot analyze id: $id1\n";
	}
	if($id2 =~ m/(chr\d):(\d+)\-(\d+)/){
		$chr2 = $1;
		my $start = $2;
		my $end = $3;
		$pos2 = int(($start + $end)/2);
	}else{
		die "cannot analyze id: $id1\n";
	}

	if($chr1 ne $chr2){
		if($chr1 eq 'chr1'){
			if($chr2 eq 'chr2'){
				$SUM_chr1_chr2 += $value;
				$INTER_chr1_chr2++;
			}elsif($chr2 eq 'chr3'){
				$SUM_chr1_chr3 += $value;
				$INTER_chr1_chr3++;
			}
		}elsif($chr1 eq 'chr2'){
			if($chr2 eq 'chr1'){
				$SUM_chr1_chr2 += $value;
				$INTER_chr1_chr2++;
			}elsif($chr2 eq 'chr3'){
				$SUM_chr2_chr3 += $value;
				$INTER_chr2_chr3++;
			}
		}elsif($chr1 eq 'chr3'){
			if($chr2 eq 'chr1'){
				$SUM_chr1_chr3 += $value;
				$INTER_chr1_chr3++;
			}elsif($chr2 eq 'chr2'){
				$SUM_chr2_chr3 += $value;
				$INTER_chr2_chr3++;
			}
		}
	}else{
		my $distance = abs($pos1 - $pos2);
		$SUM{$chr1}{$distance} += $value;
		$NUM{$chr1}{$distance}++;

		if($MAX_distance{$chr1} < $distance){
			$MAX_distance{$chr1} = $distance;
		}
	}
}
$fh_map->close();





#---------------------------------------
# output average
#---------------------------------------
my $fh_out = IO::File->new($FILE_average, 'w') or die "cannot write $FILE_average: $!";
$fh_out->print("Type\tDistance\tobserved\tCombination\tAverage\n");
foreach my $chr(qw(chr1 chr2 chr3)){
	for(my $i = 0; $i <= $MAX_distance{$chr}; $i += $Resolution){
		my $Observe = exists $SUM{$chr}{$i} ? $SUM{$chr}{$i} : 0;
		my $Combination = exists $NUM{$chr}{$i} ? $NUM{$chr}{$i} : 0;
		my $value = 0;
		if($Combination > 0){
			$value = $Observe / $Combination;
		}
		$fh_out->print("$chr\t$i\t$Observe\t$Combination\t$value\n");
	}
}
$fh_out->printf("chr1_chr2\t$MAX_distance{'chr1'}\t%d\t%d\t%f\n", $SUM_chr1_chr2, $INTER_chr1_chr2, $SUM_chr1_chr2 / $INTER_chr1_chr2);
$fh_out->printf("chr1_chr3\t$MAX_distance{'chr1'}\t%d\t%d\t%f\n", $SUM_chr1_chr3, $INTER_chr1_chr3, $SUM_chr1_chr3 / $INTER_chr1_chr3);
$fh_out->printf("chr2_chr3\t$MAX_distance{'chr1'}\t%d\t%d\t%f\n", $SUM_chr2_chr3, $INTER_chr2_chr3, $SUM_chr2_chr3 / $INTER_chr2_chr3);
$fh_out->close();


#---------------------------------------
# output parameter
# y ~ a + b*exp(-c*x) + d*exp(-e*x)
#---------------------------------------
my $getParameter = 'fit$m$getPars()';
my $Rcommand = <<"END";
R --vanilla --slave <<'EOF'
 data <- read.table("$FILE_average", header=T, sep="\t")
 attach(data)

 CalcCoeficient <- function(Ave, Dis){
	 FITTING_MAX <- max(Dis[Dis <= (max(Dis) - 1000000)])
	 if(FITTING_MAX < 2.3e6){
		FITTING_MAX <- 2.3e6
	 }
	 y <- Ave[Dis < FITTING_MAX & Dis > 0 & Ave != 0]
	 x <- Dis[Dis < FITTING_MAX & Dis > 0 & Ave != 0]
	 a0 <- min(y)
	 c0 <- log((y[1]-a0)/(y[10]-a0))/(x[10]-x[1])
	 b0 <- (y[10]-a0)/exp(-x[10]*c0)/11
	 e0 <- log((y[1]-a0-b0*exp(-c0*x[1])) / (y[3]-a0-b0*exp(-c0*x[3]))) / (x[3] - x[1])
	 d0 <- (y[1] - a0 - b0 * exp(-c0*x[1])) / exp(-e0*x[1])
	 fit <- nls(y ~ a + b*exp(-c*x) + d*exp(-e*x), start=c(a=a0, b=b0, c=c0, d=d0, e=e0))
	 a <- {$getParameter}["a"]
	 b <- {$getParameter}["b"]
	 c <- {$getParameter}["c"]
	 d <- {$getParameter}["d"]
	 e <- {$getParameter}["e"]

	 ymax <- max(Ave)
	 plot(Dis, Ave, ylim=c(0,ymax), xlab='', ylab='average frequency')
	 par(new=T)
	 new <- data.frame(x=Dis)
	 prediction <- predict(fit, new)
	 plot(Dis, prediction, ylim=c(0,ymax), type='l', col='blue', axes=F, xlab='', ylab='')
	 legend("topright",lty=1, lwd=2, col="blue", paste(sprintf("%.3f",a0), "+",
 		sprintf("%.2f",b0), "*exp(-", sprintf("%.2e",c0), "*x)+",
 		sprintf("%.2f",d0), "*exp(-", sprintf("%.2e",e0), "*x)", sep=""))
	 paste(a, b, c, d, e, sep=",")
 }

 pdf("$FILE_graph")
 par(mfrow=c(3,1), oma=c(1,1,4,1), mar=c(3,4,0,0))
 c1 <- CalcCoeficient(Average[Type=="chr1"], Distance[Type=="chr1"])
 c2 <- CalcCoeficient(Average[Type=="chr2"], Distance[Type=="chr2"])
 c3 <- CalcCoeficient(Average[Type=="chr3"], Distance[Type=="chr3"])
 gabarge <- dev.off()
 paste(c1, c2, c3, sep=":")
EOF
END

my $parameters = `$Rcommand`;
$parameters =~ s/\r?\n//;
if($parameters =~ m/"(.+)"/){
	$parameters = $1;
}

{
	my $fh_par = IO::File->new($FILE_parameter, 'w') or die "cannot write $FILE_parameter: $!";
	my @params_chr = split /:/, $parameters;
	my @chr = qw(chr1 chr2 chr3);
	my @mark = qw(a b c d e);
	for(my $i = 0; $i < 3; $i++){
		my $chr = $chr[$i];
		my @p = split /,/, $params_chr[$i];
		for(my $j =0; $j < 5; $j++){
			$fh_par->print("$chr\t$mark[$j]\t$p[$j]\n");
		}
	}
	$fh_par->printf("inter\tchr1_chr2\t%f\n", $SUM_chr1_chr2 / $INTER_chr1_chr2);
	$fh_par->printf("inter\tchr1_chr3\t%f\n", $SUM_chr1_chr3 / $INTER_chr1_chr3);
	$fh_par->printf("inter\tchr2_chr3\t%f\n", $SUM_chr2_chr3 / $INTER_chr2_chr3);
	$fh_par->close();
}








